% Figures for the paper
% 2017-2-27
% Tentative plan: We provide figures for houston
clc; clear all; close all;
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');

% predX_M1: covariates at the tract centroid is also estimated via
% hierachical model 
% M1: includes intercept and log(1+dist), coefficients are varying with Area
% M2: includes intercept and log(1+dist) and log(dist 2 coast), coefficients are varying with Area
estXmode = 3; %1=M1, 2=M2, 3=M3

% !!!!
% ONLY WORKS FOR TABLE (NOT FIGURES)
% NOW FIGURES ARE UPDATED TOO (5/25/2017)

%% Figure option
cityid = 3360; %520 for Atlanta, 1600 for Chicago, 4480 LA, 5600 NY, Houston 3360
nyear = 6;
modelind = 6; %2: linear, 6: quartic, size-cubic
islog  = 2; %2 for log(1+x)

do_fig1a = 0;
do_fig1b = 1;
do_fig1b_o = 1; %original level
do_fig1c = 1;
do_fig1c_o = 1; %original level
do_tab   = 0; %land value table

%% Load Data and preliminary stuff
load('data_all');
% load('results_full_logp1.mat');
% load('results_full_logp1_A.mat');
load('results_full_logp1_A_s3_REML.mat');
load('results_estX_REML.mat');
if estXmode == 1
    predX = estX_M1.predX;
    predM = estX_M1;
elseif estXmode == 2
    predX = estX_M2.predX;
    predM = estX_M2;
elseif estXmode == 3
    predX = estX_M3.predX;
    predM = estX_M3;
end


% other info
ncity = size(ZZ2,1);

% Spatial function
ngrid = 100;
xgrid = linspace(0,60,ngrid)';
logxgrid = log(1+xgrid);

logxgrid1 = log(1+xgrid);
logxgrid4 = [logxgrid1, logxgrid1.^2, logxgrid1.^3, logxgrid1.^4];



%% Figure 1-b: 2-d polynomial (Linear versus quartic)
if do_fig1b == 1
    % load estimation results
    workpath = pwd;
    
    % Get Sptial function using "get_sp2"
    
    % model to generate figures
    temp_M1 = est_s2; %linear
    temp_M2 = est_s6; %quartic
    
    % temp_bet = fixedEffects(temp_M2.model);
    YD = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
    
    temp_bet  = temp_M2.bet;
    temp_bet2 = temp_M2.bet2;
    temp_bet3 = temp_M2.c10;
    resid_fe = YY - XX*temp_bet - [XX(:,14).^2, XX(:,14).^2]*temp_bet2 ...
        -logDC*temp_bet3;
    
    % Actual figure
    L1 = 'Covariate-adjusted Log Price';
    L2 = 'Linear';
    L3 = 'Quartic';
    L4 = 'Integration points (tract centers)';
    xlabel_w = 'Miles from City Hall';
    ylabel_w = 'Price per Acre (1000$)';
    figtitlehead = 'Model Log(1+x)-spatial functions :';
    
    fignamehead = 'fig_2d_log_';
    % figpath = [workpath, filesep, 'figure', filesep, 'fig08_all_logp1'];
    figpath = [workpath, filesep, 'figure_draft_rev', filesep, 'fig01_ind'];
    chk_dir(figpath);
    
        
%     for cind = 1:1:ncity

        fig2 = figure;
        setmyfig(fig2, [2, 0.2, 10, 6.5])
        
%         citynum = cind;
%         cityid = xtval(cind);

            cityid_set = (1:1:pmsa_n)';
            citynum = cityid_set(xtval==cityid);
            cind = citynum;
        
        temp_tract = tract.min_dist(tract.combofips==xtval(cind),:);
        temp_YY   = YY(xtflag == xtval(cind),:);
        temp_XX   = XX(xtflag == xtval(cind),:);
        temp_grid = DD(xtflag == xtval(cind),:);
        
        temp_resid = resid_fe(xtflag==xtval(cind),:);
        
        %plot(temp_grid, temp_resid, '*', 'linewidth',2);
%         scatter(temp_grid,temp_resid,'filled', ...
%             'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('lightgrey'));
        
        scatter(temp_grid,temp_resid,'filled', ...
            'MarkerFaceColor',rgb('lightgrey'));
        
        hold on
        
        % first line
        temp_val0 = zeros(size(logxgrid1,1), nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, size(xgrid,1),1), logxgrid1];
            
            temp_val0(:,tind) = temp_C*temp_M1.post.m(:,cind);
        end
        temp_val = mean(temp_val0,2);
        %     temp_val = temp_val0;
        
        plot(xgrid,temp_val,'-*','linewidth', 2, 'Markersize', 7, 'color', rgb('grey'))
        %plot(xgrid,temp_sp1(:,cind),'-*','linewidth', 2, 'Markersize', 7, 'color', rgb('grey'))
        
        % second line
        temp_val0 = zeros(size(logxgrid4,1), nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, size(xgrid,1),1), logxgrid4];
            
            temp_val0(:,tind) = temp_C*temp_M2.post.m(:,cind);
        end
        temp_val = mean(temp_val0,2);
        %     temp_val = temp_val0;
        
        plot(xgrid,temp_val,'linewidth', 7, 'color', rgb('black'))
        %     plot(xgrid,temp_sp3(:,cind),'-*','linewidth', 4)
        %     plot(xgrid,temp_sp4(:,cind),'-o','linewidth', 4)
        
        axis tight
        %ylim([min(temp_resid), max(temp_resid)])
        ylim([9,18]);
        temp_ylim = ylim;
        plot(temp_tract,temp_ylim(1)*ones(size(temp_tract,1),1),'x','Linewidth', 2, 'Markersize',10, 'color', rgb('black'))
        hold off
        %     title([figtitlehead,name.pmsaname{name.pmsafips==cityid},', combofips = ', num2str(cityid)])

        xlim([min(temp_tract), max(temp_tract)]);
        
        xlabel(xlabel_w);
        ylabel(ylabel_w);
        
        %l = legend(L1,[L2,', ',num2str(mse1,'%4.2f')],[L3, ', ', num2str(mse2,'%4.2f')],[L4, ', ', num2str(mse3,'%4.2f')],L5);
        l = legend(L1,L2,L3,L4);
        set(l, 'fontsize', 17);
        set(gca,'fontsize', 20, 'linewidth',2);
        
        ax = gca;
        ax.YTick = log([25, 100, 800, 6400]*1000);
        ax.YTickLabel = {'25', '100', '800', '6400'};

       
        
        % save
        cd(figpath);
        saveas(gca, [fignamehead,num2str(cityid,'%04i'), '.png']);
        cd(workpath);
        
%         close
%     end
end

%% Figure 1-b: 2-d polynomial (Linear versus quartic), original level
if do_fig1b_o == 1
    % load estimation results
    workpath = pwd;
    
    % Get Sptial function using "get_sp2"
    
    % model to generate figures
    temp_M1 = est_s2; %linear
    temp_M2 = est_s6; %quartic
    
    % temp_bet = fixedEffects(temp_M2.model);
    YD = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
    
    temp_bet  = temp_M2.bet;
    temp_bet2 = temp_M2.bet2;
    temp_bet3 = temp_M2.c10;
    resid_fe = YY - XX*temp_bet - [XX(:,14).^2, XX(:,14).^2]*temp_bet2 ...
        -logDC*temp_bet3;
    
    % Actual figure
    L1 = 'Covariate-adjusted Price';
    L2 = 'Linear';
    L3 = 'Quartic';
    L4 = 'Integration points (tract centers)';
    xlabel_w = 'Miles from City Hall';
    ylabel_w = 'Price per Acre ($)';

    figtitlehead = 'Model Log(1+x)-spatial functions :';
    
    fignamehead = 'fig_2d_level_';
    % figpath = [workpath, filesep, 'figure', filesep, 'fig08_all_logp1'];
    figpath = [workpath, filesep, 'figure_draft_rev', filesep, 'fig01_ind'];
    chk_dir(figpath);
    
%     for cind = 1:1:ncity
        fig2 = figure;
        setmyfig(fig2, [2, 0.2, 10, 6.5])
        
%         citynum = cind;
%         cityid = xtval(cind);
        
            cityid_set = (1:1:pmsa_n)';
            citynum = cityid_set(xtval==cityid);
            cind = citynum;
        
        temp_tract = tract.min_dist(tract.combofips==xtval(cind),:);
        temp_YY   = YY(xtflag == xtval(cind),:);
        temp_XX   = XX(xtflag == xtval(cind),:);
        temp_grid = DD(xtflag == xtval(cind),:);
        
        temp_resid = resid_fe(xtflag==xtval(cind),:);
        
        %plot(temp_grid, temp_resid, '*', 'linewidth',2);
        scatter(temp_grid,exp(temp_resid),'filled', ...
            'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('lightgrey'));
        hold on
        
        % first line
        temp_val0 = zeros(size(logxgrid1,1), nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, size(xgrid,1),1), logxgrid1];
            
            temp_val0(:,tind) = exp(temp_C*temp_M1.post.m(:,cind));
        end
        temp_val = mean(temp_val0,2);
        %     temp_val = temp_val0;
        
        plot(xgrid,temp_val,'-*','linewidth', 2, 'Markersize', 7, 'color', rgb('grey'))
        %plot(xgrid,temp_sp1(:,cind),'-*','linewidth', 2, 'Markersize', 7, 'color', rgb('grey'))
        
        % second line
        temp_val0 = zeros(size(logxgrid4,1), nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, size(xgrid,1),1), logxgrid4];
            
            temp_val0(:,tind) = exp(temp_C*temp_M2.post.m(:,cind));
        end
        temp_val = mean(temp_val0,2);
        %     temp_val = temp_val0;
        
        plot(xgrid,temp_val,'linewidth', 7, 'color', rgb('black'))
        %     plot(xgrid,temp_sp3(:,cind),'-*','linewidth', 4)
        %     plot(xgrid,temp_sp4(:,cind),'-o','linewidth', 4)
        
        axis tight
        %ylim([min(temp_resid), max(temp_resid)])
        %     ylim([9,18]);
        temp_ylim = ylim;
        plot(temp_tract,temp_ylim(1)*ones(size(temp_tract,1),1),'x','Linewidth', 2, 'Markersize',10, 'color', rgb('black'))
        hold off
        %     title([figtitlehead,name.pmsaname{name.pmsafips==cityid},', combofips = ', num2str(cityid)])
        
        xlim([min(temp_tract), max(temp_tract)]);
        
        xlabel(xlabel_w);
        
        %l = legend(L1,[L2,', ',num2str(mse1,'%4.2f')],[L3, ', ', num2str(mse2,'%4.2f')],[L4, ', ', num2str(mse3,'%4.2f')],L5);
        l = legend(L1,L2,L3,L4);
        set(l, 'fontsize', 15);
        set(gca,'fontsize', 15, 'linewidth',2);
        
        ylim([0, 2000000])
        
        
        % save
        cd(figpath);
        saveas(gca, [fignamehead,num2str(cityid,'%04i'), '.png']);
        cd(workpath);
        
%         close
%     end
end

%% Figure 1-c: 3-d  land values, log scale
if do_fig1c == 1
    m_XX = mean(XX);
    m_XX2 = mean([XX(:,14).^2, XX(:,14).^3]);
    XX2 = [XX(:,14).^2, XX(:,14).^3];
    
    eval(['temp_est = est_s', num2str(modelind),';']);
    temp_sig2 = temp_est.model.MSE; %residual error variance
    
    figtitlehead = 'Land Value at the tract center :';
    
    workpath = pwd;
    fignamehead = 'fig_3d_log_';
    figpath = [workpath, filesep, 'figure_draft_rev', filesep, 'fig01_ind_estX', num2str(estXmode)];
    chk_dir(figpath);
    
%     for j=1:1:length(xtval)
        fig3 = figure;
        setmyfig(fig3, [2, 0.2, 10, 6.5])
        
            cityid_set = (1:1:pmsa_n)';
            citynum = cityid_set(xtval==cityid);
            j = citynum;
        
%         citynum = j;
%         cityid = xtval(j);
        
        fig = figure;
        setmyfig(fig, [1.7, 0.2, 10, 6.5],3)
        
        % temporary variable (for transaction)
        temp_XX  = XX(xtflag == xtval(j),:); %covariates in this city
        temp_XX2 = XX2(xtflag == xtval(j),:);
        temp_DD  = DD(xtflag == xtval(j),:);
        temp_ZZ  = ZZ(xtflag == xtval(j),:);
        
        % temporary variable (for tract)
        cityid = xtval(j);
        temp_sel  = tract.combofips==cityid;
        temp_grid = tract.min_dist(temp_sel);
        temp_lon = tract.lon(temp_sel);
        temp_lat = tract.lat(temp_sel);
        temp_d2c = tract.d2c(temp_sel);
        temp_n = size(temp_grid,1);
        
        temp_log_grid = log(1+temp_grid);
        temp_grid4 = [temp_log_grid, temp_log_grid.^2, temp_log_grid.^3, temp_log_grid.^4];
        
        
        % contructing distance dependent covariates for prediction
        % predX: we estimate and predict outside of this m-file.
        temp_predX1 = predX(temp_sel,1:14);
        temp_predX2 = predX(temp_sel,15:16);
        
        % land value at traact center in time t
        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*log(temp_d2c) + temp_est.post.sig2/2; % value does not depend on time
        
        % value that depends on time
        temp_val = zeros(temp_n, nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, temp_n,1), temp_grid4];
            temp_val(:,tind) = temp_val_nv + temp_C*temp_est.post.m(:,j) ...
                + 1/2*diag(temp_C*temp_est.post.V(:,:,j)*temp_C');
        end
        %temp_val = exp(temp_val);
        temp_val = mean(temp_val,2); % log scale
        
        
        %scatter3(temp_lon, temp_lat,(temp_val), 'filled', 'r')
%         scatter3(temp_lon, temp_lat,(temp_val), 'filled', ...
%             'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('orangered'));
        
        scatter3(temp_lon, temp_lat,(temp_val), 'filled', ...
            'MarkerFaceColor',rgb('orangered'));
        
        hold on
        
        % 3-d surface on the grid
        temp_x = temp_lon;
        temp_y = temp_lat;
        %temp_z = temp_val;
        temp_z = (temp_val);
        
        tri = delaunay(temp_x,temp_y);
        % plot(temp_x,temp_y,'.')
        [r,c] = size(tri);
        % disp(r)
        h = trisurf(tri, temp_x, temp_y, temp_z);
        axis vis3d
        % axis off
        % l = light('Position',[-50 -15 29]);
        % set(gca,'CameraPosition',[208 -50 7687])
        lighting phong
        shading interp
        colormap gray
        % shading flat
        % colorbar EastOutside
        hold off
        set(gca, 'fontsize', 15, 'linewidth', 3);
        xlabel('Longitude');
        ylabel('Latitude');
        zlabel('Land values');
        %     title([figtitlehead,name.pmsaname{name.pmsafips==cityid},', combofips = ', num2str(cityid)])
                
        ax = gca;
        ax.ZTick = log([100, 200, 400, 800]);
        ax.ZTickLabel = {'100', '200', '400', '800'};
        
        % save
        cd(figpath);
        saveas(gca, [fignamehead,num2str(cityid,'%04i'), '.png']);
        cd(workpath);
        close all
        
%     end
end

%% Figure 1-c: 3-d  land values
if do_fig1c_o == 1
    m_XX = mean(XX);
    m_XX2 = mean([XX(:,14).^2, XX(:,14).^3]);
    XX2 = [XX(:,14).^2, XX(:,14).^3];
    
    
    eval(['temp_est = est_s', num2str(modelind),';']);
    temp_sig2 = temp_est.model.MSE; %residual error variance
    
    figtitlehead = 'Land Value at the tract center :';
    
    workpath = pwd;
    fignamehead = 'fig_3d_level_';
    figpath = [workpath, filesep, 'figure_draft_rev', filesep, 'fig01_ind_estX', num2str(estXmode)];
    chk_dir(figpath);
    
%     for j=1:1:length(xtval)
        fig3 = figure;
        setmyfig(fig3, [2, 0.2, 10, 6.5])
        
            cityid_set = (1:1:pmsa_n)';
            citynum = cityid_set(xtval==cityid);
            j = citynum;
        
%         citynum = j;
%         cityid = xtval(j);
        
        fig = figure;
        setmyfig(fig, [1.7, 0.2, 10, 6.5],3)
        
        % temporary variable (for transaction)
        temp_XX  = XX(xtflag == xtval(j),:); %covariates in this city
        temp_XX2 = XX2(xtflag == xtval(j),:);
        temp_DD  = DD(xtflag == xtval(j),:);
        temp_ZZ  = ZZ(xtflag == xtval(j),:);
        
        % temporary variable (for tract)
        cityid = xtval(j);
        temp_sel  = tract.combofips==cityid;
        temp_grid = tract.min_dist(temp_sel);
        temp_lon = tract.lon(temp_sel);
        temp_lat = tract.lat(temp_sel);
        temp_d2c = tract.d2c(temp_sel);
        temp_n = size(temp_grid,1);
        
        temp_log_grid = log(temp_grid);
        temp_grid4 = [temp_log_grid, temp_log_grid.^2, temp_log_grid.^3, temp_log_grid.^4];
        
        
% contructing distance dependent covariates for prediction
        % predX: we estimate and predict outside of this m-file.
        temp_predX1 = predX(temp_sel,1:14);
        temp_predX2 = predX(temp_sel,15:16);
        
        % land value at traact center in time t
        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*log(temp_d2c) + temp_est.post.sig2/2; % value does not depend on time
        
        % value that depends on time
        temp_val = zeros(temp_n, nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, temp_n,1), temp_grid4];
            temp_val(:,tind) = temp_val_nv + temp_C*temp_est.post.m(:,j) ...
                + 1/2*diag(temp_C*temp_est.post.V(:,:,j)*temp_C');
        end
        %temp_val = exp(temp_val);
        temp_val = mean(exp(temp_val),2)/1000; % level scale
        
        
        %scatter3(temp_lon, temp_lat,(temp_val), 'filled', 'r')
%         scatter3(temp_lon, temp_lat,(temp_val), 'filled', ...
%             'MarkerFaceAlpha',4/8,'MarkerFaceColor',rgb('orangered'));
        
        scatter3(temp_lon, temp_lat,(temp_val), 'filled', ...
            'MarkerFaceColor',rgb('orangered'));
        
        hold on
        
        % 3-d surface on the grid
        temp_x = temp_lon;
        temp_y = temp_lat;
        %temp_z = temp_val;
        temp_z = (temp_val);
        
        tri = delaunay(temp_x,temp_y);
        % plot(temp_x,temp_y,'.')
        [r,c] = size(tri);
        % disp(r)
        h = trisurf(tri, temp_x, temp_y, temp_z);
        axis vis3d
        % axis off
        % l = light('Position',[-50 -15 29]);
        % set(gca,'CameraPosition',[208 -50 7687])
        lighting phong
        shading interp
        colormap gray
        % shading flat
        % colorbar EastOutside
        hold off
        set(gca, 'fontsize', 15, 'linewidth', 3);
        xlabel('Longitude '); %(^{\circ})
        ylabel('Latitude'); %(^{\circ})
        zlabel('Price per Acre (1000$)');
        %     title([figtitlehead,name.pmsaname{name.pmsafips==cityid},', combofips = ', num2str(cityid)])
        
%         zlim([0, 6400]);
        ax = gca;
        ax.ZTick = ([25, 100, 200, 400, 800, 1600]);
        ax.ZTickLabel = {'25', '100', '200', '400', '800', '1600'};        
        
        % save
        cd(figpath);
        saveas(gca, [fignamehead,num2str(cityid,'%04i'), '.png']);
        cd(workpath);
%         close all
        
%     end
end

